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Abstract. We apply the Gradient-Newton-Galerkin-Algorithm (GNGA) of Neuberger & Swift 
to find solutions to a semilinear elliptic Dirichlet problem on the region whose boundary is the 
Koch snowflake. In a recent paper, we described an accurate and efficient method for generating a 
basis of eigenfunctions of the Laplacian on this region. In that work, we used the symmetry of the 
snowflake region to analyze and post-process the basis, rendering it suitable for input to the GNGA. 
The GNGA uses Newton's method on the eigenfunction expansion coefficients to find solutions to 
the semilinear problem. This article introduces the bifurcation digraph, an extension of the lattice 
of isotropy subgroups. For our example, the bifurcation digraph shows the 23 possible symmetry 
types of solutions to the PDE and the 59 generic symmetry-breaking bifurcations among these 
symmetry types. Our numerical code uses continuation methods, and follows branches created 
at symmetry-breaking bifurcations, so the human user does not need to supply initial guesses for 
Newton's method. Starting from the known trivial solution, the code automatically finds at least one 
solution with each of the symmetry types that we predict can exist. Such computationally intensive 
investigations necessitated the writing of automated branch following code, whereby symmetry 
information was used to reduce the number of computations per GNGA execution and to make 
. intelligent branch following decisions at bifurcation points. 
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1. Introduction. 

We seek numerical solutions to the semilinear elliptic boundary value problem 

Au + f x (u) = in Q 

in', (I) u = on an, 

O ■ 

where A is the Laplacian operator, $7 C R 2 is the region whose boundary dQ is the Koch snowflake, 
u : 0, — > R is the unknown function, and fx : R — > R is a one-parameter family of odd functions. For 
convenience, we refer to as the Koch snowflake region. This article is one of the first to consider 
a nonlinear PDE on a region with fractal boundary. In this paper, we choose the nonlinearity to 
be 
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(2) f x (u) = \u + ti 
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and treat A 6 R as the bifurcation parameter. When the parameter is fixed, we will sometimes use 



/ in place of f\. Using this convention, note that A = /'(0). 

This paper exploits the hexagonal symmetry of the Koch snowflake region, and the fact that / 
is odd. Our nonlinear code would work with any region with hexagonal symmetry and any odd 
'superlinear' function / (see [3]), and with minor modification for other classes of nonlinearities 
as well. We chose to work with odd / primarily because of the rich symmetry structure. The 
explicit shape of ft represents a considerable technological challenge for the computation of the 
eigenfunctions |16} 127] . which are required as input to the nonlinear code. 

It is well known that the eigenvalues of the Laplacian under this boundary condition satisfy 

(3) < Ai < A 2 < A 3 < >oo, 
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and that the corresponding eigenfunctions {V'jljeN can be chosen to be an orthogonal basis for 
the Sobolev space H = Hq(Q) = W ' 2 (Q), and an orthonormal basis for the larger Hilbert space 
L 2 = L 2 (Q). The inner products are 

(u,v)h = / Vu-Vvdx and (u,v)2 = / uvdx, 
Jn Jn 

respectively (see [TJ [151 E!)- Theorem 8.37 and subsequent remarks in [9] imply that the 
eigenfunctions are in C°°(J2). In [T7J, properties of the gradients of eigenfunctions near boundary 
points are explored in light of the lack of regularity of d£l. 

Using the Gradient-Newton-Galerkin- Algorithm (GNGA, see [26]) we seek approximate solutions 
u = 2 i=i a j^j t° © by applying Newton's method to the eigenfunction expansion coefficients of 
the gradient V J{u) of a nonlinear functional J whose critical points are the desired solutions. The 
definition of J, the required variational equations, a description of the GNGA, and a brief history 
of the problem are the subject of Section [2j 

The GNGA requires as input a basis spanning a sufficiently large but finite dimensional subspace 
Bm = spanj^i, . . . , V'm}, corresponding to the first M eigenvalues As described in [27] . 

a grid Gm of N carefully placed points is used to approximate the eigenfunctions. These are the 
same grid points used for the numerical integrations required by Newton's method. Section [3] briefly 
describes the process we use for generating the eigenfunctions. 

Section U] concerns the effects of symmetry on automated branch following. The symmetry the- 
ory for linear operators found in [27] is summarized and then the extensions required for nonlinear 
operators are described. Symmetry-breaking bifurcations are analyzed in a way that allows an 
automated system to follow the branches created at the bifurcations. As we develop the theory, we 
present specific examples applying the general theory to equation ([1]) on the snowflake region. In 
particular, we find that there are 23 different symmetry types of solutions to (pQ), and 59 generic 
symmetry-breaking bifurcations. The symmetry types and bifurcations among them are summa- 
rized in a bifurcation digraph, which generalizes the well-known lattice of isotropy subgroups (see 
|10|). As far as we know, the bifurcation digraph is a new way to organize the information about 
the symmetry-breaking bifurcations. 

Section [5] describes how understanding the symmetry allows remarkable increases in the efficiency 
of the GNGA. Section [6] describes the automated branch following. We use repeated executions of 
the GNGA or a slightly modified algorithm (parameter-modified GNGA) to follow solution branches 
of ([U E]). The GNGA uses Newton's method, which is known to work well if it has a good initial 
approximation. The main shortcoming of Newton's method is that is works poorly without a good 
initial approximation. We avoid this problem by starting with the trivial solution (u = 0). The 
symmetry-breaking bifurcations of the trivial solution are found by the algorithm and the primary 
branches are started. The program follows the branches by continuation methods, and then follows 
the new branches created at symmetry-breaking bifurcations. To follow an existing branch, we vary 
A slightly between executions. To start new solution branches created at bifurcation points, we 
treat A as a variable and fix one of the null eigenfunctions of the Hessian evaluated at the bifurcation 
point. The symmetry analysis tells which null eigenfunction to use. In this way solutions with all 
23 symmetry types are found automatically, starting from u = 0, without having to guess any 
approximations for Newton's method. 

In our experiments, many bifurcation diagrams were generated by applying the techniques men- 
tioned above. A selection of these diagrams are provided in Section [7J along with contour plots 
of solutions to ([T]) corresponding to each of the 23 symmetry types predicted to exist. We include 
evidence of the convergence of our algorithm as the number of modes M and grid points N increase. 

Many extensions to our work are possible, including enforcing different boundary conditions on 
the same region, solving similar semilinear equations on other fractal regions, and applying the 
methodology to partial difference equations (PdE) on graphs [25] . Section [8] discusses some of 
these possible extensions. In particular, we are in the process of re-writing the suite of programs. 
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We plan to be able to solve larger problems using a parallel environment. We will be able to solve 
problems with larger symmetry groups by automating the extensive group theoretic calculations. 
This concluding section also has a discussion of the convergence of the GNGA. 

2. GNGA. 

We now present the variational machinery for studying ([I]) and follow with a brief description of 
the general GNGA. Section [6] contains more details of the implementation of the algorithm for our 
specific problem. Let F\(u) = J Q f\(s)ds for all u G R define the primitive of f\. We then define 
the action functional J : R x H — > R by 

(4) J(X,u) = [ {±\Vu\ 2 - F x (u)} dx. 

Jn 

We will sometimes use J : H — > R to denote J(A, •). The class of nonlinearities / found in 
[33 O [23| |2"5] imply that J is well defined and of class C 2 on H . The choice ([5]) we make in 
this paper belongs to that class. Critical points of J are by definition weak solutions of ([JJ (see 
for example [UIMllS]), and clearly classical solutions are critical points. The usual "bootstrap" 
argument of repeatedly applying Theorem 8.10 of [U] can be used in our case. Specifically, Hq is 
embedded in L q for all q > 2 when the space diminsion n is 2, regardless of the regularity of d£l (due 
to the zero Dirichlet boundary condition, see PQ). Hence u G H k implies f(u) G H k as well. As a 
result, if u is a critical point then u G C°°(0) n C(f2), hence a classical solution. If one considered 
boundary conditions, space dimensions, and nonlinear terms other than the choices made in this 
paper, it could happen that critical points would be weak not classicial solutions. Regardless, our 
approximations lie in Bm C C°°. Here, the existence proofs for positive, negative, and sign-changing 
exactly once solutions from [U [2H] immediately give at least 3 nontrivial (classical) solutions for 
our specific superlinear boundary value problem; appealing to symmetry implies the existence of 
even more solutions (see for example [25J). 

The choice of H for the domain is crucial to the analysis of the PDE (see [IJ [24] , and references 
therein), as well as for understanding the theoretical basis of effective steepest descent algorithms 
(see [7[I22[[23]. for example). We will work in the coefficient space R M = Bm- The coefficient vector 
of u G Bm is the vector a G R M satisfying u = Y^j=i a j' l l J j ■ Using the corresponding eigenvalues 
([3]) and integrating by parts, the quantities of interest are 

(5) gj = J'(u)(^-) = / {Vu ■ Vvpj - f(u) ipj} = ajXj - / f(u) ipj, and 

Jn Jn 

(6) h jk = J"(u)(ipj,^ k ) = / {Vtjjj ■ Vvpk - f'(u) ipj i) k } = Xj5 jk - / f(u) ipj ip k , 

Jn Jn 

where 5jk is the Kronecker delta function. Note that there is no need for numerical differentiation 
when forming gradient and Hessian coefficient vectors and matrices in implementing Algorithm 12.11 
this information is encoded in the eigenfunctions. 

The vector g G R M and the M xM matrix h represent suitable projections of the I? gradient and 
Hessian of J, restricted to the subspace Bm, where all such quantities are defined. For example, 
for u = J2jL\ %"0i> v = YlfLi bjipj, and w = ^j=i Cj^ji we have: 

M 

Pb m V2«7(m) = ^2 dj^ji J'( u )( v ) = 9 ' b, and J"(u)(v, w) = Kb ■ c = b ■ he. 
j=i 

We can identify g with the approximation -Pb m V2 J(u) of V2 J{u) = An + f(u), which is defined 
for u G Bm- The solution x to the M-dimensional linear system h\ = 9 is then identified with the 
(suitably projected) search direction (D\ J(u)) _1 V2«/(m), which is not only defined for u G Bm, but 
is there equal to {D 2 H J(n))~ 1 Vn J(u). We use the least squares solution of h\ = g. In practice, 
the algorithm works even near bifurcation points where the Hessian is not invertible. 
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The heart of our code is Newton's method in the space of eigenfunction coefficients: 

Algorithm 2.1. (GNGA) 

(1) Choose initial coefficients a = {a,j}jL l , and set u = ^2,ajipj. 

(2) Loop 

(a) Calculate the gradient vector g = {J' (u)(ipj)}^L 1 from equation (J5]) . 

(b) Calculate the Hessian matrix h = {J"(u)(ipj, ipk)'}j I fi=i ^ rom equation ([6]). 

(c) Exit loop if 1 1^| | is sufficiently small. 

(d) Solve hx = 9 for the Newton search direction x £ ^ M • 

(e) Replace a a — x an d update u = ^ a>ji/)j ■ 

(3) Calculate sig(/i) and J for the approximate solution. 

If Newton's method converges then we expect that u approximates a solution to the PDE (fT]), 
provided M is sufficiently large and the eigenfunctions and numerical integrations are sufficiently 
accurate. See Section 

Our estimate for the Morse index (MI) of the critical point of J is the signature of h, denoted 
sig(/i), which is defined as the number of negative eigenvalues of h. This measures the number of 
linearly independent directions away from u in which J decreases quadratically. 

The basic Algorithm 12.11 is modified to take advantage of the symmetry of our problem. The M 
integrations required in step (a) and the M{M + l)/2 integrations in step (b) are reduced to fewer 
integrations if the initial guess has nontrivial symmetry. 

We often use a "parameter- modified" version of the GNGA (pmGNGA). In this modification, A 
is treated as an unknown variable and one of the M coefficients is fixed. Along a given branch, 
symmetry generally forces many coefficients to be zero. When a bifurcation point is located by 
observing a change in MI, we can predict the symmetry of the bifurcating branches using the 
symmetry of the null eigenfunctions of the Hessian. By forcing a small nonzero component in the 
direction of a null eigenfunction (orthogonal to the old branch's smaller invariant subspace), we can 
assure that the pmGNGA will not converge to a solution lying on the old branch. Another benefit 
of the pmGNGA is that it can handle a curve bifurcating to the right as well as one bifurcating 
to the left. In our system, the branches that bifurcate to the right have saddle node bifurcations 
where they turn around and go to the left. The pmGNGA can follow such branches while the 
normal GNGA cannot. 

The implementation of pmGNGA is not difficult. The M equations are still 

gi = J'(u)(ipi) = 0, 

but the M unknowns are 

a = (ai, . . . , a,k-i, A, a,k+i, • • • , (im), 
and the value of one coefficient, a&, is fixed. Consequently, we replace the Hessian matrix h with a 
new matrix h where the fe-th column is set to dgi/dX = —ai: 

niJ ~\ ~a t \ij = k ■ 
The search direction x is the solution to the system hx = g- The pmGNGA step is 

a <~ a - X, 

and then u and A are updated. After Newton's method converges, the A;-th column of the original 
is calculated and the MI of the solution, sig(/i), is computed. 
We conclude this section with a very brief history of the analytical and numerical aspects of 
the research into (pQ) given our type of nonlinearity /. Our introduction to this general subject 
was [4], where a sign-changing existence result was proven. This theorem is extended in [5]; we 
indicate briefly in Section [7] where this so-called CCN solution can be found on our bifurcation 
diagrams. The article [7] was our first success in using symmetry to find higher MI solutions. The 
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Figure 1. The Koch snowflake region Q with the grids G13 and G133 at levels 
£ = 2 and 3, respectively. A generic grid point (which is not on any line of reflection 
symmetry) is indicated in the larger grid. 

GNGA was developed in [26 , wherein a much more detailed description of the variational structure 
and numerical implementation can be found. The first implementation of the GNGA for regions 
where the eigenfunctions are not known in closed form is in |12j . where the region is a Bunimovich 
stadium. The article [24J provides a historical overview of the authors' experimental results using 
variants of the Mountain Pass Algorithm (MPA, MMPA, HLA) and the GNGA, as well as recent 
analytical results and a list of open problems; the references found therein are extensive. 



3. The Basis of Eigenfunctions. 

In [27], we describe theoretical and computational results that lead to the generation of a basis 
of eigenfunctions solving 

(7) An + An = in 0, u = on dQ. 

That paper details the grid technique and symmetry analysis that accompanied the effort; we briefly 
summarize those results in this section. 

The Koch snowflake is a well-known fractal, with Hausdorff dimension log 3 4. Following Lapidus, 

Neuberger, Renka, and Griffith |16j . we take our snowflake to be inscribed in a circle of radius ^ 
centered about the origin. We use a triangular grid Gn of N points to approximate the snowflake 
region. Then, we identify u : Gn — > M. with u 6 M. N , that is, 

(8) u(xi) = Ui 

at grid points X{ £ Gn- Our paper [27] differs from [16] in that we use a different placement of the 
grid points and a different method of enforcing the boundary condition, resulting in more accurate 
eigenvalue estimates with fewer points. Figured] depicts the levels 2 and 3 grids in the family of grids 
used in [27] to compute eigenfunctions; we used the first M eigenfunctions computed at levels 4, 5, 
and 6 in our nonlinear experiments. The number of grid points at level I is N = (9^ — 4^) /5, and the 
spacing between grid points is h = 2/3^. We computed the eigenvalues and eigenfunctions for ([7]) 
using ARPACK and this approximation to the Laplacian with zero-Dirichlet boundary conditions: 

(9) -A«(i) ~ ((12 — number of neighbors) u(x) — ^^{neighboring values of u}^ . 

The ARPACK is based upon an algorithmic variant of the Arnoldi process called the Implicitly 
Restarted Arnoldi Method (see [19]) and is ideally suited for finding the eigen-pairs of the large 
sparse matrices associated with the discretization of the Laplacian. 
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4. Symmetry: The Lattice of Isotropy Subgroups and The Bifurcation Digraph. 

This section describes equivariant bifurcation theory as it applies to the branching of solutions 
to equation ([1]), see [H [THl [HI HHJ- We are able to describe the expected symmetry types of 
solutions to (HJ, as traditionally arranged in a lattice of isotropy subgroups. We introduce the 
bifurcation digraph, a refinement of the lattice, which shows every possible generic bifurcation from 
one symmetry type to another as a directed edge which is labeled with information about the 
bifurcation. The bifurcation digraph is of interest in its own right and summarizes the essential 
information required by our automated branch following code. In this project, GAP (Groups, 
Algorithms, and Programming, see [8]) was used solely to verify the symmetry analysis we did by 
hand. In our continuing projects GAP is a useful tool since it can perform the tedious calculations 
and write the information in a format that can be read by the branch following code. Matthews 
|21j has used GAP to do similar calculations. We apply this methodology to the snowflake domain 
being considered in this paper. The analysis shows that solutions fall into 23 symmetry types, and 
that there are 59 types of generic symmetry breaking bifurcations. 

Group Actions and the Lattice of Isotropy Subgroups. Let T be a finite group and V be a 
real vector space. A representation of T is a homomorphism a : T — > GL(V). Where convenient, 
we identify GL(V) with the set of invertible matrices with real coefficients. Every representation 
a corresponds to a unique group action of T on V by the rule 7 • v := 0(7) (v) for all 7 G T and 
v G V. We will usually use the action rather than the representation. The group orbit of v is 
r • v = {7 • v I 7 G r}. 

Example 4.1. Let 

D 6 := (p,a\ p 6 = a 2 = 1, pa = ap 5 } 

be the dihedral group with 12 elements. It is convenient to define r = p 3 a. It follows that 
err = tct = p 3 . The group H>q is the symmetry of a regular hexagon, and of the Koch Snowflake 
region Q. The standard ID>6 action on the plane is given by 

p ■ (x, v) = (h x + 1? y> -ir* + b) 

(!0) cr ■ (x, y) = (-x, y) 

t- (x,y) = (x,-y). 

In this action, p is a rotation by 60°, a is a reflection across the y-axis, and r is a reflection across 
the x-axis. These group actions are depicted in Figure [131 near the end of the paper. 

We will denote subgroups of ©6 by listing the generators. While any given subgroup of ©6 
can be defined using only p and a, we find it geometrically descriptive to use r in certain cases. 
For example, we prefer (p 2 ,r) to the equivalent (p 2 ,po~). In order to make relationships among 
subgroups intuitive, we often include r when its membership is implied by the other generators 
(see for example Figure [2]) . 

The standard B)q group action (|10p is not the only action we consider. For a function u G L 2 (£l) 
and group element 7 G B)q, we define (7 • u)(x) = u{^~ 1 ■ x). In this paper, a vector u defined by 
Ui = u(xi), for a given grid Gn = {xi}^L 1 , is a discrete approximation of a function on Q. The ~D>q 
group action on u G is a permutation of the components: (j-u)i = u(7 _1 - xi). Given a function 
u G L 2 (Q) or WL N , the group orbit ID>6 • u consists of functions obtained from u by a reflection or 
rotation. 

Example 4.2. The group x Z2, where Z2 = {1,-1}, acts on L 2 (Q.) in a natural way. For all 

(7, z) G IH>6 x ^2, define 

(7, z) ■ u = z(j • u). 

We will denote (7, 1) G Dg x Z2 by 7 and (7, —1) G ©6 x ^2 by —7. With this natural notation 
(—7) ■ u = — (7 • u), which we call simply —7 • u. 
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Let us recall some facts about group actions, following [61 \TU[ lllj . The isotropy subgroup or 
stabilizer of v G V in F is 

Stab(u, r) := {7 e T | 7 • -u = v}. 

The isotropy subgroup measures the symmetry of v, and is sometimes called the little group of v, or 
r„. If the group T is understood, we may simply write Stab(u) in place of Stab(i>, T). The stabilizer 
of a subset W C V in T is Stab(PV, T) := {7 G T | 7 • W = W}. This must be distinguished from 
the point stabilizer of a subset 

P Stab(Ty, r) := {7 G T | 7 • v = v for all v G W} = p|{Stab(v, r) | v G TV}. 

Another commonly used notation is r^y for the stabilizer and IV^n for the point stabilizer. Note 
that pStab(W, T) is always normal in Stab(W, T), and the effective symmetry group acting on W is 
Stab(W, T)/ pStab(VK, T) , which acts faithfully on W . 

If £ is a subgroup of V then the fixed point subspace ofT, inV is 

Fix(S, V) := {v G V | 7 • v = v for all 7 e S}. 

Another notation for the fixed point subspace is Vs. We write Fix(S) when V is understood. 

An isotropy subgroup of the V action on V is the stabilizer of some point v G V. For some group 
actions, not every subgroup of V is an isotropy subgroup. 

Example 4.3. Consider the B>q action on the plane M? described in equation (jlOi It is well-known 
that (p) is not an isotropy subgroup of this action. 

Now consider the Hq action on the function space L 2 (Q). We give a standard argument that 
every subgroup of H>q is an isotropy subgroup. Start with a function u* that is zero everywhere 
except for a small region, and suppose that the region is distinct from each of its nontrivial images 
under the B>q action. Then for any subgroup S < U)q, the average of the function u* over S, defined 
as 

(11) PZ(U*) = T^J2^- U * 

has isotropy subgroup S. Therefore every subgroup of the B>q action on L 2 (ft) is an isotropy 
subgroup. The average over the group is an example of a Haar operator, and Ps : V — > Fix(S, V) 
is an orthogonal projection operator [36] . 

Similarly, every subgroup of Hq is an isotropy subgroup of the Bg action on 1R , the space of 
functions on the grid Gat, provided I > 3. This follows from averaging the function that is 1 at a 
generic lattice point, and elsewhere. Recall that a generic point is one whose isotropy subgroup 
is trivial. Figure Q] shows that the level two grid G13 does not have a generic point, while the level 
three grid G133 does. Thus, the space of functions on G133 has the same isotropy subgroups as 
L 2 (il), but a much smaller space has this same property. Start with any generic point x\ G f2. 
Then U)q acts on the space of functions on the 12 points Og • x±. This B>q action on M 12 has the 
same structure of isotropy subgroups as the Dg action on L 2 (0), and is the Bg action used in 
our GAP calculations. The corresponding 12-dimensional representation is the well-known regular 
representation of Bg (see [29 ] [31 ] 134"]). 

The symmetry of functions is described by two related concepts. A function q : V — > R is 
T-invariant if 5(7 • v) = q(v) for all 7 G T and all v G V. Similarly, an operator T : V — > V is 
T-equivariant if T(7 • v) = 7 ■ T(v) for all 7 G T and all v G V. 

Example 4.4. The energy functional J defined in equation @ is Bg x ^-invariant. The nonlinear 
PDE dl]) can be written as (A + f)(u) = 0, where A + / is a Bg x Z2-equivariant operator. (There 
are subtleties concerning the domain and range of A. See [6l [7] for a careful treatment of the 
function spaces.) In particular, A + / is Bg-equivariant because the snowflake region Q has Bg 
symmetry, and (A + f)(—u) = —(A + f)(u), since / is odd. As a consequence, if u is a solution to 
(H|, then so is every element in its group orbit (Bg X Z2) • u. 
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The isotropy subgroups and fixed point subspaces are important because of the following simple 
yet powerful results. See [6} fTOl [TT] . 

Proposition 4.5. Suppose V acts linearly on V , T : V —> V is T-equivariant and £ is an isotropy 
subgroup ofT. 

(a) Ifve Fix(E) then T(v) G Fix(E). Thus, T| Fix(E) : Fix(£) -»• Fix(£) is defined. 

(b) Stab(Fix(£)) = iVr(E), the normalizer ofT, in T, and pStab(Fix(£)) = £. 

(c) T| Fix ( S ) is Nr(S)-equivariant. 

(d) T|Fi x (£) ?s iVr(E) /H-equivariant, anrfiVr(E)/E acts faithfully on Fix(E). 

If E is a subgroup of T, the normalizer of E in Y is defined to be A r r(E) := {7 G T | 7E = E7}, 
which is the largest subgroup of V for which E is a normal subgroup. The presence of the normalizer 
in Proposition 14.5( b) is interesting, since the normalizer is a property of the abstract groups, and 
is independent of the group action. 

Example 4.6. As a consequence of Proposition l3~5l we can solve the PDE ([I]), written as (A+/)(u) = 
0, by restricting u to functions in Fix(E, L 2 (S7)). This leads to a simpler problem since the function 
space Fix(£,L 2 (Q)) is simpler than L 2 (Q). An example of this is in Costa, Ding, and Neuberger 
[7]. The techniques of that paper, applied to our problem, would find sign-changing solutions with 
Morse index 2 within the space Fix(ID>6, L 2 (Q)). This space consists of all functions which are 
unchanged under all of the rotations and reflections of the snowflake region. 

Proposition (j4.5f) also applies to the GNGA, since the Newton's method iteration mapping is 
H)q x Z2-equivariant. If the initial guess is in a particular fixed point subspace, all the iterates will 
be in that fixed point subspace. This fact can be used to speed numerical calculations, as described 
in Section 

Two subgroups £i,£2 of T are conjugate (£1 ~ £2) if £1 = 7£27~ 1 for some 7 G T. The 
symmetry type of v G V for the T action is the conjugacy class of Stab(t> , T) . Note that Stab(7 • v) = 
7Stab(u)7 _1 . Thus, every element of a group orbit V ■ v has the same symmetry type. 

Let S = {Si} denote the set of all symmetry types of a T action on V. The set S has a natural 
partial order, with Si < Sj if there exits £» G Si and Ej G Sj such that £» < E,-. The partially 
ordered set (5, <) is called the lattice of isotropy subgroups of the T action on V [10]. The diagram 
of the lattice of isotropy subgroups is a directed graph with vertices Si and arrows Si <— Sj if, and 
only if, Si < Sj and there is no symmetry type between Si and Sj. 

Example 4.7. The symmetry type of a solution u to our PDE (pQ) for the Dq x 7*2 action is the 
conjugacy class of Stab(u, Dq x Z2); we refer to this as the symmetry type of n, without reference 
to Dg x Z2. The discussion of Dq acting on L 2 (VL) in Example 14.31 can easily be extended to x Z2 
acting on L 2 (Cl). Note that if — 1 G £ < B>q x Z2, then the average of any function over £ is u = 0. 
Therefore the only isotropy subgroup of Oq x 7*2 which contains — 1 is Dq x Z2 itself. On the other 
hand, the argument in Example 14.31 shows that any subgroup of ©6 x ^2 which does not contain 
— 1 is an isotropy subgroup. Therefore, £ < Dg x Z2 is an isotropy subgroup of this group action 
if and only if £ = D 6 x Z 2 or -1 $ £. 

This result allowed us to compute the isotropy subgroups by hand. We verified our calculations 
using GAP. There are exactly 23 conjugacy classes of isotropy subgroups for the ~Bq x Z2 action 
on L 2 (f2), shown in condensed form in Figure [2j Thus, a solution to the PDE ([TJ has one of 23 
different symmetry types. 

Irreducible Representations and the Isotypic Decomposition. In order to understand the 
symmetry-breaking bifurcations we need to first understand irreducible representations and the 
isotypic decomposition of a group action. The information about the irreducible representations is 
summarized in character tables [291 EH E21 El] • For our purposes, irreducible representations over 
the field M are required, see [6l 1101 ITT]. The irreducible representations ofT are homomorphisms 
from r to the space of dj x dj real matrices: 7 1— > 0^(7), such that no proper subspace of M dj is 
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x Z, 
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r 6 = 


{-a, -r) 


r 7 = 


(°, -r) 


r 8 = 





to 




to 




r 15 = 


{<?) 


r 16 = 


(r) 


Tir = 


<-r> 


r 18 = 


(-a) 



Ti = (p,cr,r) = ] 
T 2 = (p, -cr, -t) 
T3 = {-p,a,-r) 
T 4 = {-p,-a,r) 



r 9 = 




Tio 






= (p 2 ,-r) 


Tl2 







Figure 2. The condensed diagram of the isotropy lattice (see [10]) for the H)q x Z2 
action on L 2 (f2). The vertices of this diagram are the symmetry types (equivalence 
classes of isotropy subgroups). We follow the convention [6], [lOl [TT] that one element 
T.; of each symmetry type Si = [Ti] is listed. The representatives Tj have the property 
that Ti < Tj iff S, 



< Sj. Contour plots of solutions to PDE ((TJ) with each of the 



23 symmetry types are given in Figures [13] and Q3] The diagram of the isotropy 
lattice is condensed as in |32j. The small numbers on the edges tell the number of 
connections emanating from each symmetry type in a box. A missing small number 
means 1. For example, the two arrows representing fr 2 i] < [^13] and [F 21 ] < [r 14 ] 
in the full diagram are collapsed to a single arrow in the condensed diagram. For 
Tq through T^, the r generator is redundant since r = p^a, but its presence makes 



the subgroups manifest. For example, T 2 = (p, 
generators make it clear that {—cr, — r) < {p, —a, 



-t) 



-t) = {p,—a), but the three 
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invariant under a^(j) for all 7 G T. The dimension of the irreducible representation a^' is dj. We 
call W C V a T -invariant subspace of V if T • C W. An irreducible subspace of V is an invariant 
subspace with no proper invariant subspaces. Every irreducible subspace of the V action on V 
corresponds to a unique (up to similarity) irreducible representation of T. The dimension of the 
irreducible subspace is the same as the dimension of the corresponding irreducible representation. 

For each irreducible representation a^' of V, the isotypic component of V for the T action, 

(i) 

denoted by Vp , is defined to be the direct sum of all of the irreducible subspaces corresponding 
to the fixed a^' [BJ [TUl EJ E7j. The isotypic decomposition of V is then 

(12) V = ®V®. 

j 

Some of the isotypic components might be the single point at the origin. These can be left out 
of the isotypic decomposition. A description of the isotypic components in terms of projection 
operators is given in |27j . 

For any group T, we denote the trivial representation by That is a^(j) = 1 for all 7 G T. 
Thus, if T is an isotropy subgroup of a To action on V, then 

y r (1) = Fix(r,F). 

Example 4.8. Let us consider the = (p, o~, r) action on L 2 (il). We need to consider the six 
irreducible representations of H)q, which are listed in [27], to find the isotypic decomposition of 
L 2 (Q). Since these isotypic components are central to our problem, we drop the U)q an d define 



(13) 



vJ£ , j = 1, 2, . . . , 6 as follows: 



yd) 


= {u 




p- 


u = 


u, a ■ u = u, t • u = u} 




y( 2 ) 


= {u 


£ L 2 (n) 


p- 


u = 


u, a ■ u = —u, t ■ u = - 


-«} 


y( 3 ) 


= {u 


£ L 2 (Q) 


p- 


u = 


—u, a ■ u = u, t ■ u = - 


-«} 


y( 4 ) 


= {u 


e L 2 (fl) 


p- 


u = 


—u, a • u = —u, t ■ u = 


u} 




= {u 


G L 2 (n) 


P 3 


• u = 


= u, u + p 2 ■ u + p 4 ■ u = 


0} 


y( 6 ) 


= {u 


g L 2 (n) 


P 3 


• u = 


- —u, u + p 2 ■ u + p 4 ■ u 


= 0} 



Example 4.9. The isotypic decomposition of Ti3 = (p) 
sentation theory. The irreducible representations of Z, 



o 



(P) 



t/3)J- 



representations of 
(14) 



1 for j = 1, 2, . . . , 6. Over the field 
Zq are given by 

1 



= 7*Q illustrates some features of real repre- 
over C are all one-dimensional. They are 
., however, the one-dimensional irreducible 



a 



(l) 



(P) 



a 



(2) 



(P) 



and the two-dimensional irreducible representations of 
given by 



-1, 



up to similarity transformations, are 



(15) 



a 



(3) 



(p) 




a 



('!) 



(P) 



l 
2 

V3 
2 

Note that a^(p) is matrix for a rotation by 120° and a^(p) is a 60° rotation matrix. 

An irreducible representation over R is called absolutely irreducible if it is also irreducible over C. 
For example, all of the irreducible representations of H)q listed in (27J are absolutely irreducible, as 
are the one-dimensional irreducible representations of in equation (|14p . On the other hand, the 
two-dimensional irreducible representations of 7L§ in equation (|15p are not absolutely irreducible. 

The four isotypic components of the (p) action on L 2 (J7) are 

,(i) 



V, 



(i>) 



{u G L 2 {VL) j p ■ u = u} = V {1) V {2) 
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V 



,(2) 




If we had used the complex irreducible representations, some of the corresponding isotypic compo- 
nents would contain complex-valued functions. It is more natural to use real irreducible representa- 
tions, and consider only real-valued functions. The price we pay is that most of the representation 
theory found in books, and built into GAP, is done for complex irreducible representations. 

The isotypic decomposition for each of the 23 isotropy subgroups, Tj, of Dq x Z2 can be written 
as a direct sum of some subset of the eight spaces for j = 1,...,4, and v} j) and V 2 U) for 

j = 5,6 defined in (113p and |27| . The C++ program can easily check if a function is in any of the 

(7) 

isotypic components of Bm for each of the Tj, i = 0, 1, . . . , 22, actions. 

Symmetry-Breaking Bifurcations. The fact that there are 23 possible symmetry types of 
solutions to the PDE (pQ) begs the question, do solutions with each of these symmetry types exist? 
Clearly the trivial solution u = 0, with symmetry type So, exits. Our procedure for finding 
approximate solutions with each of these symmetry types is to start with the trivial solution and 
recursively follow solution branches created at symmetry-breaking bifurcations. 

Let us start by abstracting the PDE defined by ([T]), which depends on the real parameter A. Let 
V be an inner product space and J:Mxy->lbea family of To— invariant functions that depends 
on a parameter A. That is, J(A, 7 • u) = J(A, u) for all 7 G Tq and u G V. It is understood that Tq 
is the largest known group for which J is invariant; of course J is also invariant under any subgroup 
of Tq. We will use T, or Tj, to refer to an isotropy subgroup of the "full" group Tq. Consider the 
steady-state bifurcation problem g(X,u) = 0, where g(X,u) = VJ(X,u). Throughout this paper, 
the gradient V acts on the u component. The solutions to g{X, u) = are critical points of J, so 
we use the terms "solution" and "critical point" interchangeably. Note that g:M.xV—>Vis& 
family of To— equivariant gradient operators on V. That is, g(X, 7 • u) = 7 • g(X,u). For our PDE, 
Tq =IJiq x 7j2- Ln the numerical implementation, V = = Bm and g is defined in ([5j). 

We define a branch of solutions to be a connected component of {(X,u) elx L 2 (f2) | g(X,u) = 
0, Stab(u) = T}, where T is called the isotropy subgroup, or symmetry, of the branch. A branch of 
solutions B\ has a symmetry-breaking bifurcation at the bifurcation point (X*,u*) G B\ if a branch 
of solutions, B2, with a different symmetry, has (X*,u*) as a limit point but (X*,u*) £ Bi- We say 
that branch B2 is created at this bifurcation, and often refer to B\ as the mother branch and B2 
as the daughter branch. The symmetry of the daughter branch is always a proper subgroup of the 
symmetry of the mother branch. That is, the daughter has less symmetry than the mother. 

The main tool for finding bifurcation points is the Hessian of the energy functional, h. If (X* ,u*) 
is a bifurcation point, then h(X*,u*) is not invertible, since otherwise the implicit function theorem 
would guarantee the existence of a unique local solution branch. The Morse index (MI) of a critical 
point (A, it) is defined to be the number of negative eigenvalues of h(X,u) = D 2 J(X,u), provided 
no eigenvalue is 0. The Hessian is symmetric, so all of its eigenvalues are real. The MI on a branch 
of solutions typically changes at a bifurcation point. 

Example 4.10. The trivial solution to (HJ [2]) is u = 0, and the trivial branch is {(A, 0) | A G M}. 
Since h(X,0)(v) = Av + Xv, the bifurcation points of the trivial branch are (Aj,0), where A,,i G N, 
are the eigenvalues (J3J) - If Aj < A < Aj+i, then the MI of the trivial solution (A,0) is i. The i-th 
primary branch is created at the bifurcation point (Aj, 0) on the trivial branch. In cases with double 
eigenvalues there are two branches created at the same point in our problem. For example, the 
second and third primary branches are created at A2 = A3. Near (Aj,0), the solutions on the i-th 
primary branch are approximately some constant times the i-th eigenfunction of the Laplacian, fy. 

We define a degenerate critical point, or a degenerate solution, to be a point (A*, it*) which 
satisfies g(X*,u*) = and det h(X*, u*) = 0. Thus, every bifurcation point is a degenerate critical 
point. Some degenerate critical points are not bifurcation points. For example, when a branch 
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folds over and is not monotonic in A, the fold point is degenerate, but is not a bifurcation point as 
we have defined it. (Note that we avoid the term "saddle- node bifurcation" since there is really no 
bifurcation.) 

Let us develop some notation to talk about bifurcations. Suppose that (A*, it*) is an isolated 
degenerate critical point of a To-equivariant system g(X,u) = 0. Let T = Stab(u*,ro), and define 
L := h(\*,u*). Note that T, not Tq, is important as far as the bifurcation of (A*, it*) is concerned. 
Let E be the null space of the T-equivariant operator L. We call E the center eigenspace. Let T' 
be the point stabilizer of E. The definitions are repeated here for reference: 

(16) r := Stab(n*,r ), L := h(X*,u*), E := N(L), V := P Stab(£,r). 

If e E E, then L(e) = by definition. For any 7 S T, 7 • e £ E since the T-equivariance of L 
implies that L(j • e) = 7 • L(e) = 0. Hence, 

Stab(£,r) = r. 

Note that St&b(E,T)/ pStab(£, T) = T/V acts faithfully on E. In the usual case where (A*,u*) is 
a bifurcation point, not just a degenerate critical point, we say that T/V is the symmetry group of 
the bifurcation, or that (A*, it*) undergoes a bifurcation with T/V symmetry. 

In the notation of ([To]) . L sends each of the isotypic components Vp to itself [271 ED GS]- 
Barring "accidental degeneracy," the center eigenspace E is a T-irreducible subspace. Thus, is 
typically a subspace of exactly one isotypic component \ and dim.{E) is the dimension dj of 
the corresponding corresponding irreducible representation, Furthermore, the point stabilizer 
of E is the kernel of a^' and can be computed without knowing E. In summary, at a generic 
bifurcation point there is some irreducible representation aft' of T such that: 

E is T-irreducible, E QVf?\ dhn(E) = AMI = dj, V = {7 e T | a {j) (~f) = I}. 

Accidental degeneracy is discussed in (27J EU [3l] . We did not encounter any accidental degeneracy 
in our numerical investigation of ([HE]), so we will not discuss it further here. 

We finally have the background to describe the bifurcations which occur in equivariant systems. 
The goal is to predict what solutions will be created at each of the symmetry breaking bifurcations, 
and know what vectors in E to use to start these branches using the pmGNGA. While such a 
prediction is impossible for some complicated groups, we can determine how to follow all of the 
bifurcating branches in the system ([11 [2]) . We follow the treatment and notation of \10\ [Tl] . At a 
symmetry-breaking bifurcation we can translate (A*, it*) to the origin, and we could, in principle, 
do an equivariant Liapunov-Schmidt reduction or center manifold reduction to obtain reduced 
bifurcation equations g : R x E — > E where 5(0, 0) = 0, Dg(0, 0) = 0, and g is T := Stab(u*)- 
equivariant. It is important to realize that we do not actually need to perform the Liapunov-Schmidt 
reduction. 

The most powerful tool for understanding symmetry breaking bifurcations is the Equivariant 
Branching Lemma. Recall that absolutely irreducible representations were defined in Example 14.91 
See [U [TU1 flT] for a thorough discussion of the Equivariant Branching Lemma, including further 
references. 

Theorem 4.11. Equivariant Branching Lemma (EBL) Suppose T acts absolutely irreducibly 
on the space E, and let g : R x E — > E be T -equivariant. Assume that T acts nontrivially, so 
g(X,0) = 0. Since T acts absolutely irreducibly, Dg(X,0) = c(X)Id for some function c : R — > R, 
where Id is the identity matrix of size d = dim(.E'). Assume that c(0) = and c'(0) 7^ 0. Let X 
be an isotropy subgroup of the T action on E with dimFix(£,.E) = 1. Then there are at least two 
solution branches of g{X,u) = with isotropy subgroup S created at (0,0). 

The EBL, combined with Liapunov-Schmidt theory, implies that there are at least two solution 
branches of the full problem g(A, u) = with isotropy subgroup S created at the bifurcation point 
(A*, it*). We call these newly created branches EBL branches since their existence can be predicted 
by the EBL. Other branches created at a bifurcation are called non-EBL branches. 
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(p,a,r) (p,tr,r) (p,<t,t) (p,<t,t) (p,<t,t) {p,v,t) 



(p) (p 2 ,a) <p 2 ,r> (a,r) (a) 



(r) 



(P 3 



(1) 



Figure 3. Diagrams of the six isotropy lattices for the actions of Oq = (p, a, r) on 
each of the six isotypic components V^O of the B6 action on L 2 (0). This describes 
the six possibilities (barring accidental degeneracy) for the Bg action on the center 
eigenspace E at a degenerate critical point. 

Following [6tllO|ITT]. we define a maximal isotropy subgroup of a V action on V to be an isotropy 
subgroup E/r with the property that if is an isotropy subgroup such that E < 0, then = E 
or = r. In other words, a maximal isotropy subgroup is a maximal proper isotropy subgroup. If 
dim(Fix(E, E)) = 1, then E is a maximal isotropy subgroup of the T action on E. The converse, 
however, is not true. 

In gradient systems, for example the PDE ([1]), more can be said. If E is any maximal isotropy 
subgroup of the V action on E, then there is typically a solution branch created at the bifurcation 
with isotropy subgroup E. If dimFix(E,.E) > 2, the branch created is an example of a non-EBL 
branch. See [30J for a discussion of bifurcations in gradient systems. 

By Proposition 14. 5| the effective symmetry group of g, restricted to Fix(E,.E), is iVr(E)/E. This 
effective symmetry group determines how solutions with symmetry E bifurcate. 

Example 4.12. Consider a degenerate critical point with isotropy subgroup Ti = ©6 = (p,a,r). 
Barring accidental degeneracy, the center eigenspace £ is a subspace of one of the 6 isotypic 
components of the Bg action on L 2 (Q) described in Example 14.81 Figure [3] shows the lattice of 
isotropy subgroups for B)q acting on each of these 6 isotypic components V^'. These 6 cases can 
be distinguished by determining which isotypic component an arbitrary eigenfunction in E belongs 
to. We shall go through each of these six cases, and describe the resulting bifurcation. Recall that 
r = T\ = 3q for each of these six cases, and T' = pStab(£',r). 



E C V {1) = 


> r' 


= r 1 = 


= {p,v,t), 


dimE = 


1, 


r/r' = 


(i) 


E C V {2) = 


> r' 


= r 13 


= </»>, 


dimE = 


1, 


r/r'^ 


z 2 


E C V {3) = 


> r' 


= r 9 = 


= (p 2 ,v), 


dimE = 


1, 


r/r' ^ 


z 2 


E C = 


> r' 


= r 10 


= (p 2 ,r), 


dimE = 


1, 


r/r' ^ 


z 2 


E C T/ {5) = 


> r' 


= r 19 


= (P 3 ), 


dimE = 


2, 


r/r' ^ 


B 3 


E C V {6) = 


> r' 


= T22 


= (1), 


dimE = 


2. 


r/r' ^ 


B 6 



The first case, E C = Fix(ri, L 2 (il)), does not lead to a symmetry-breaking bifurcation. The 
Bg action on E is trivial, so the EBL does not apply. The degenerate critical point (n*,A*) is 
typically a fold point (or saddle-node), not a bifurcation point. In the neighborhood of the fold 
point there is only one solution branch, with isotropy subgroup V±, and the branch lies to one side 
of A = A* or the other. 

The next three cases, with T /V = Z 2 symmetry, are called pitchfork bifurcations. Clearly, the 
only maximal isotropy subgroup is T' in each case, and the EBL applies. The effective symmetry 
group acting on E is Z 2 , so there are two conjugate solution branches created at the bifurcation. 
In the branch following code we follow one of these branches using the pmGNGA starting with any 
eigenvector e £ E. 



11 



JOHN M. NEUBERGER, NANDOR SIEBEN, AND JAMES W. SWIFT 



(p) (p) (P) (P) 



(P 2 ) (P 3 ) (1) 

Figure 4. The diagrams of the four isotropy lattices for the actions of Ti3 = (p) on 
each of the four isotypic components V$ of the Ti3 action on L 2 (S7). This describes 
the four possibilities (barring accidental degeneracy) for the Ti3 action on the center 
eigenspace E at a degenerate critical point. 

The next case, with E C is a bifurcation with D3 symmetry. The maximal isotropy 

subgroup = (a, r) satisfies 

dimFix(r 5 ,£) = 1, and N ri {T B )/T 5 = (1). 

Our branch following code uses a projection operator to find an eigenvector e £ E with Stab(e, T±) = 
T5. The pmGNGA using this eigenvector e will follow one of the solution branches created at the 
bifurcation, and the pmGNGA using the negative eigenvector — e will find a branch that is not 
conjugate to the first. Bifurcations with B3 symmetry are typically transcritical, and two D3-orbits 
of branches are created at the bifurcation [TOj E] . 

The last case, with E C y( 6 ), is a bifurcation with Dg symmetry. There are two maximal 
symmetry types, the conjugacy classes of and r^. A calculation shows that 

dimFix(r 15 ,£) = dimFix(r 16 ,£) = 1, and N Tl (T 15 )/T 15 = N Tl (T 16 )/T 16 = Z 2 . 

To follow one branch from each of the group orbits of solution branches created at this bifurcation, 
it suffices to use the pmGNGA twice, with the eigenvectors ei,e 2 6 E, where Stab(ei,ri) = Ti5 
and Stab(e2,ri) = T\q. It is well-known that these EBL-branches are typically the only branches 
created at a bifurcation with ©6 symmetry [10} [11] . 

Example 4.13. Consider a degenerate critical point with isotropy subgroup Ti3 = (p) = Zq. Barring 
accidental degeneracy, the center eigenspace E is a subspace of one of the 4 isotypic components 
V$ defined in Example 14.91 Figure H] shows the lattice of isotropy subgroups for Ti3 acting on 

each of these 4 isotypic components. Recall that T = Tx3 = (p) for each of these cases, and the 
minimal isotropy subgroup is V = pStab(£', T). We shall go through each of the four cases, and 
describe the resulting bifurcation: 

v (i) eV (2) ^ r' = r 13 = {p), dim£ = i, r/r'^(i) 

v (3) eV (4) ^ y' = T 2X = {p 2 ), di mJ e = i, r/r'^z 2 

r' = r 19 = (p 3 ), dims = 2, r/r'^z 3 

r' = r 22 = (i), dims = 2, r/r'^Zg. 

The first two cases are analogous to the first two cases in Example 14. 121 When T/V = (1) there is a 
fold point, but no symmetry breaking bifurcation. There is a pitchfork bifurcation when T/T' = Z 2 . 
The next two cases are interesting because Ti3 does not act absolutely irreducibly on E, and the 
EBL does not apply. In both cases T' is a maximal isotropy subgroup. 

In the third case, where E C V$ = V® , every eigenfunction in the 2-dimensional E has 
isotropy subgroup Tig. Since we have a gradient system, we know that solution branches with 
isotropy subgroup Tig are created at this bifurcation with Z3 symmetry. The bifurcation is well- 
understood, and it looks like a bifurcation with B3 symmetry, except that the "angle" of the 
bifurcating solutions in the E plane is arbitrary. This means that trial and error is needed, in 
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general, to find eigenfunctions in E for which the pmGNGA will converge. If a branch is found for 
a starting eigenfunction e, then the eigenfunction — e is used to find the other solution branch. 

In the fourth case, where E C V^j = every eigenfunction in E has the same isotropy 

subgroup: T22 = (!)■ Gradient bifurcations with "Lq symmetry look like bifurcations with B@ 
symmetry, except that the angle in the E plane is arbitrary. Again, trial and error is needed to 
find starting eigenfunctions for which the pmGNGA converges. 

The Bifurcation Digraph. A calculation similar to those summarized in Examples 14.121 and 
14.131 was done for each of the isotropy subgroups of the H>q x Z2 action on L 2 (f2). The calculations 

were done by hand, and verified with GAP. There are 59 generic symmetry-breaking bifurcations, 

(7) 

one for each isotypic component on which acts nontrivially. The amount of information is 
overwhelming, so we display the essential results in what we call a bifurcation digraph. 

Definition 4.14. The bifurcation digraph of the Tq action on a real vector space V is a directed 
graph with labelled arrows. The vertices are the symmetry types (equivalence classes of isotropy 
subgroups). Given £ < T, two isotropy subgroups of the Tq action on V, we draw an arrow from 
[r] to [E] iff £ is a maximal isotropy subgroup of the T action on some isotypic component 

of V. Each arrow has the label T/T', where V is the kernel of the T action on . Furthermore, 
each arrow is either solid, dashed or dotted. The arrow is 

solid if dimFix(£,£) = 1 and JV r (E)/E = Z 2 , 

dashed if dimFix(£,£) = 1 and iVr(£)/E = (1), and 

dotted if dimFix(£,£) > 2, 

(i) 

where E is any irreducible subspace contained in Vp . 

Note that if dimFix(£,.E) = 1, then iVr(£)/£ is either Z2 or (1), since these are the only linear 
group actions on E = R 1 . Thus, the three arrow types (solid, dashed, and dotted) exhaust all 
possibilities. 

Theorem 4.15. For a given Tq action on V, every arrow in the diagram of the isotropy lattice is 
an arrow in the bifurcation digraph. 

Proof. Suppose [r] —> [£] is an arrow in the diagram of the isotropy lattice. Then some £* G [E] is a 
maximal isotropy subgroup of the T action on V. Choose u* G V such that Stab(u*, T) = £*. Such 
a u* exists since E* is an isotropy subgroup. Now consider the isotypic decomposition {V^}jeJ 
of V. We can write u* = YljeJ u ^i where G are uniquely determined. Let 7 be any 
element of E*. Then 7 • u* = Sjej 7 " u<y ^ = u * ■ Since each of the components Vp is T-invariant, 
7 • U U) = u (j) for each j G J. Thus E* < Stab(u (j) , T) for each j G J. Either Stab(n (i) ,r) = T or 
Stab(w^, r) = E*, since E* is a maximal isotropy subgroup of the V action on V. If Stab(it^,r) = 

r for all j G J, then Stab(n*,r) = T. But Stab (it*, T) 7^ T, so Stab(tt^\r) = E* for some j G J, 

(i) 

and E* is a maximal isotropy subgroup of the T action on this component V-f of V. Therefore the 
bifurcation digraph has an arrow from [r] to [E*] = [E]. □ 

Theorem 14.151 says that the bifurcation digraph is an extension of the diagram of the isotropy 
lattice. The bifurcation digraph has more arrows, in general. As with the lattice of isotropy 
subgroups, we usually draw a single element T of the equivalence class [r] for each vertex of the 
bifurcation digraph. 

An arrow from T to E in the bifurcation digraph indicates that a To-equivariant gradient system 
g(X, u) = can have a generic symmetry-breaking bifurcation where a mother branch with isotropy 
subgroup r creates a daughter branch with isotropy subgroup E. The symmetry group of the 
bifurcation is T/V, and the center eigenspace at the bifurcation point is the T-irreducible space 
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E. The information encoded in the label and arrow type is used by the heuristics of our branch- 
following algorithm. A solid arrow indicates that every e in the one-dimensional space Fix(Y,,E) 
satisfies 7 • e = — e for some 7 G S. Thus, there is typically a pitchfork bifurcation in the space 
Fix(£, E). A dashed arrow indicates that 7 • e = e for all e G Fix(£,£) and 7 G X. Thus, 
the daughter branches bifurcating in the directions e and — e are not conjugate. A dotted arrow 
indicates that the EBL does not apply to this bifurcation. As mentioned above, branching of 
solutions corresponding to a dotted arrow is generic in gradient systems |30[ ITU] . 

A condensed bifurcation digraph for the Dg x Z2 action on L 2 (Q) is shown in Figure [5j The 
calculations for the directed edges coming from T\ and Ti3 are described in examples 14. 121 ane 14. 13[ 
respectively. The digraph has 65 directed edges, but there are only 5 possibilities for the symmetry 
group of the bifurcation: Y /V = Z2, Z3, Zq, B3, or Mq. The symmetry-breaking bifurcation with 
each of these symmetries is well understood [10| lllj. and each is described briefly in Example 14.121 
or 14.131 This digraph is of great help in writing an automated code for branch following. 

In our problem the label T /V and arrow type are sufficient to characterize the bifurcation 
completely. For more complicated groups, the label may need to contain more information about 
the action of T on E. For example the label T/V = §4 would be ambiguous, since §4 has two 
faithful irreducible representations with different lattices of isotropy subgroups. 

5. Symmetry and Computational Efficiency. 

Several modifications of the GNGA f)2. 1 1) take advantage of symmetry to speed up the calcula- 
tions. The symmetry forces many of the components of the gradient and Hessian to be zero. We 
identified these zero components and avoided doing the time-consuming numerical integrations to 
compute them. At the start of the C++ program, the isotropy subgroup, Tj, of the initial guess is 
computed. Recall that there are M modes in the Galerkin space Bm, so dim(Sj\,f) = M. Define 
Mi := dim(Fix(r,, Bm))- We chose the representatives Tj within each conjugacy class so that 
Fix(Fj, Bm) is a coordinate subspace of Bm- Thus, M — Mi components of the gradient g(X, u) are 
zero if Fix(tt) = Tj. The numerical integrations in ([5]) are done only for the Mi potentially nonzero 
components of g. Similarly, Mj(Mj + l)/2 rather than M(M + l)/2 numerical integrations are 
needed to compute the part of the Hessian matrix h needed by the GNGA algorithm: The numer- 
ical integrations in © are done only if tpj and V'A: are both in Fix(Fj, Bm)- The system h\ = 9 for 
the Newton step x reduces to a system of Mi equations in Mj unknowns. After Newton's method 
converges to a solution, the full Hessian needs to be calculated in order to compute the MI. Here, 
too, we can take advantage of the symmetry: Since h is Tj -equivariant, hjk = if ipj and ipk are 

(i) 

in different isotypic components Vp. of Bm- 

As an example, consider the execution time for approximating a solution with T\ symmetry 
using M = 300 modes and a level i = 5 grid on a 1GHz PC. Our C++ code uses only Mi = 30 
modes, and takes about 1.5 seconds per Newton step, compared to 44 seconds when the symmetry 
speedup is not implemented. 

6. Automated Branch Following. 

The branch following code is a complex collection of about a dozen Perl scripts, Mathematica and 
Gnuplot scripts, and a C++ program. These programs write and call each other fully automatically 
and communicate through output files, pipes and command line arguments. A complete bifurcation 
diagram can be produced by a single call to the main Perl script. 

Two choices for the function of u plotted vs. A are shown in Figure [6j In most bifurcation 
diagrams we plot approximate solutions u evaluated at a generic point (2/27,4^3/27) (the big dot 
in Figure [TJ) versus the parameter A; other choices for the vertical axis such as J(u) or ||it||oo lead to 
less visible separation of branches. Two conjugate solutions can have different values at the generic 
point, but since our program follows only one branch in each group orbit this does not cause a 
problem. 
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Figure 5. The bifurcation digraph for the H>q x Z2 action on L 2 (fi) extends the 
diagram of the isotropy lattice. The digraph shown is condensed as in Figure [2j 
The arrows indicate generic symmetry breaking bifurcations. The Morse index of 
the mother branch changes by 1 at bifurcations with Z2 symmetry, and it changes 
by 2 at all other bifurcations shown here. 



The C++ program implements the GNGA algorithm. Its input is a vector of coefficients a £ M M 
for an initial guess in Newton's method, an interval for A, a stepsize for A and several other 
parameters such as the grid level. It finds solutions on a single branch of the bifurcation diagram. 
Every solution is written as a single line in an output file. This line contains all the information 
about the solution, and can be used to write an input file for a subsequent call to the same C++ 
program. 

The C++ program finds one branch (referred to as the main branch) and a short segment of each 
of the daughter branches created at bifurcations of the main branch. The coefficients approximating 
the first solution on the branch are supplied to the C++ program. Newton's method is used to find 
this first solution, then A is incremented and the next solution is found. The program attempts to 
follow the main branch all the way to the final A, usually 0. Heuristics are used to double or halve 
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FIGURE 6. Bifurcation diagrams of the sixth primary branch (which bifurcates 
from showing \\uW2 and ii(2/27, 4\/3/27) as a function of A. Since \\uW2 is a 
©6 x ^-invariant function of u, each group orbit of solution branches is shown as 
one curve on the left. The disadvantage of plotting \ \uW2 is that the curves in many 
bifurcation diagrams are not well separated. The point (2/27, 4\/3/27) is not on 
any of the reflection axes of the snowflake region. There are 2 primary branches 
with symmetry Si, four secondary branches with symmetry Sg, and four secondary 
branches with symmetry Sio. Our choice for the bifurcation diagrams in this paper 
combines the advantages of both views: u(2/27, 4\/3/27) is plotted as a function of 
A for exactly one branch (the solid lines) from each group orbit. Unless indicated 
otherwise, all figures were produced with level I = 5 and M = 300 modes. 

the A stepsize when needed, keeping the stepsize in the interval from the initial stepsize (input to 
the C++ program) to 1/32 of the initial stepsize. For example, the stepsize is halved if Newton's 
method does not converge, if it converges to a solution with the wrong symmetry, or if more than 
one bifurcation is detected in one A step. 

The Morse index is computed at each A value on the main branch. When the MI changes a 
subroutine is called to handle the bifurcation before the main branch is continued. If the MI 
changes from mi to 777,2, we define m = max{ 7774, 777,2}. Then the bifurcation point is approximated 
by using the secant method to set the m-th eigenvalue of the Hessian h{u) to zero as a function of 
A. The GNGA is needed at each step of the secant method to compute u = u(X). We find that the 
GNGA works well even though we are approximating a solution for which the Hessian is singular. 

After the bifurcation point is approximated, a short segment of each bifurcating branch is com- 
puted and one output file is written for each branch, using Algorithm 16.11 If the Equivariant 
Branching Lemma (EBL) holds, then we know exactly which critical eigenvector to use for each 
branch. 

Algorithm 6.1. (f ollow_branch) 

(1) Input: bifurcation point (A, a) , one critical eigenvector e £ R , 

and stepsize AA < 0. Output: A file is written for one daughter branch. 

(2) Write (A, a) to output file. Set £ = 0.1. Set A b = A. 

(3) Compute index k so that |efc| > \et\ for all i 6 {1, ...,M}. 

(4) Repeat until A;, — A < AA, or t < 0.1/32 or some maximum number of points have 
been written to the file. 

(a) Do the pmGNGA with initial guess (A, a + te), fixing coefficient k. 

(b) If Newton's method converges replace (A, a) by the solution found and 
write this point to the file, else t<—t/2. 
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Figure 7. A partial bifurcation diagram of the 14-th primary branch showing a Bgj 
a D3 and several Z2 bifurcations. At the Bg bifurcation, 12 branches in two different 
group orbits are born. In accordance with Figure [H only two branches are followed 
and shown on this bifurcation diagram. An animation showing the followed branch 
with symmetry type S15 is shown in s3sl5.gif, and an animation of the followed 
branch with symmetry type Si 7 is in s3sl7s7.gif. Note that this branch with Si 7 
symmetry "dies" at a bifurcation with 7Li symmetry, showing that we cannot always 
make a consistent distinction between secondary and tertiary branches. At the O3 
bifurcation, 6 branches in two different group orbits are born. As before, only two 
branches are followed. An animation showing the "upper" branch with symmetry 
type S7, through the bifurcation point and continuing to the "lower" branch with 
symmetry type S7 is shown in s7s3s7.gif. For clarity, the branches bifurcating 
from 3 of the Z2 bifurcations are not shown. The numbers next to a branch indicate 
the MI of the solution. The MI changes by 2 at a square, and by 1 at a circle. 

Note that the pmGNGA can follow a branch that bifurcates to the right or the left. Those 
that bifurcate to the right usually turn over in a saddle-node "bifurcation" that does not offer any 
difficulty for the pmGNGA. Figures [7] and [8] show several examples of bifurcations. 

The EBL does not hold at bifurcations with Z3 and "Lq symmetry in our problem, since the 2- 
dimensional center eigenspace does not have a 1-dimensional subspace with more symmetry. Figure 
[8] shows one of the few bifurcations with Z3 symmetry that we observed. By good fortune, the 
branches with symmetry type S19 were successfully followed using the same eigenvectors one would 
choose for a bifurcation with B3 symmetry. A better method for following bifurcating solutions that 
are not predicted by the EBL would be to use the pmGNGA with random (normalized) eigenvectors 
in E repeatedly until it appears that all equivalence classes of solutions have been found. 

The branch following code is called recursively by a main Perl script. Initially, the C++ program 
follows the trivial branch on a given A range. This results in an output file for the trivial branch 
and another output file for each bifurcating primary branch. Then the short parts of the primary 
branches are followed with more calls to the C++ program. Any bifurcating branch results in 
a new output file, and the Perl script makes another call to the C++ program to continue that 
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Figure 8. The B3 bifurcation of the 13-th primary branch is on the left. This is 
the only observed B3 bifurcation that is not transcritical. An animation of the upper 
branch with symmetry type S5, through the bifurcation point and continuing with 
the lower branch is shown in s5sls5.gif. A Z3 bifurcation of a daughter of the 
24-th primary branch is shown on the right. The branches created at this bifurcation 
are not described by the EBL. An animation of the branches with symmetry type 
5*i9 is shown in sl9sl3sl9.gif. 



branch. The main Perl script's most important job is book keeping. It saves the output files with 
distinct names, and calls the branch following code to continue each of the new branches. The 
process stops when all the branches are fully followed within the given A range. 

In this way, a complete bifurcation diagram is produced by a single invocation of the main Perl 
script. There is no need to guess initial conditions for input to Newton's method, since the trivial 
solution is known exactly (a = 0) and all the other solutions are followed automatically. 

The main Perl script calls several other smaller scripts. For example, there is a script which 
extracts solutions from output files and feeds them to the branch following code as input. Another 
script creates Gnuplot scripts on the fly to generate bifurcation diagrams. Perl scripts are used to 
automatically number and store the output files and create human readable reports about them. 



7. Numerical Results. 

Our goal was to find solutions to ([TJ [2]) at A = with each of the 23 symmetry types. The 
24-th primary branch is the first one with symmetry type S2, so we followed the first 24 primary 
branches. With level £ = 5 and M = 300 modes, which gave our most accurate results, this found 
solutions with all symmetry types except S\\ and S14. We then searched the first 100 primary 
branches, only following solutions with symmetry above Sn and Si 4 on the bifurcation digraph 
(Figure [5j) In this way we found solutions with all 23 symmetry types. The bifurcation diagrams 
which lead to these solutions are shown in Figures [9H121 We chose one solution at A = with 
each symmetry type by taking the one descended from the lowest primary branch. These choices 
are indicated by dots in Figures [9 H121 and the corresponding contour diagrams of the solutions are 
shown in Figures [131 and [Til The contour diagrams use white for u > and black for u < 0, and 
gray indicates u = 0. Equally spaced contours are drawn along with dots for local extrema. Details 
about the technique for generating these contour diagrams are found in [27] . 

We ran our experiments using a range of modes and levels in order to observe convergence and 
qualitative stability of the implementation of our algorithm. At level £ = 5 we have computed 300 
eigenfunctions so M < 300 is possible. At level £ = 6 we computed only 100 eigenfunctions. Due 
to our limited computational resources, using more than 100 modes on level 6 was not practical. 
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Figure 9. The complete bifurcation diagram for the first six primary branches 
bifurcating from the trivial branch. The second branch, with symmetry £7, contains 
the CCN solution. The dots at A = in Figures [9HT21 correspond to solutions 
depicted in Figures [13] and Q31 We used the level 5 grid with 300 modes in creating 
all bifurcation diagrams. In Figure [TBI convergence data for the solution of symmetry 
type S10 at A = is provided. 



As an indication of the convergence, consider the bifurcation diagram in Figure The diagram 
looks qualitatively the same for any choice of I and M that we used. The position of the bifurcation 
point creating the S\q solution (near A = 30) changes slightly, according to this table: 





£ = 4 


I = 5 


£ = 6 


M = 


100 


35.3931 


34.9814 


34.9252 


M = 


200 


32.1131 


32.2964 




M = 


300 




32.0518 





The level 5 and 6 approximations with M = 100 modes are very close, but increasing the mode 
number has more of an effect. This indicates that the results with (£,M) = (5,300) are more 
accurate than those with (6, 100). Figure US shows how u(2/27, 4^/3/27) varies with mode number 
and i for the solution with Sio symmetry at A = shown in Figures [9] and [T3j The horizontal 
segments of the graphs correspond to the addition of modes with zero coefficients for this solution. 
Based on this and other similar convergence results, we chose to use level 5 with 300 modes in most 
of our numerical experiments. 



8. Conclusions. 

We are currently working on a more general program for recursive branch following in symmetric 
systems. Starting with any graph, the analog to Equation [JJ is the Partial difference Equation 
(PdE) Lu + f(u) = [25], where L is the well-known discrete Laplacian on that graph and u is 
a real-valued function on the vertices. Discretizing a PDE as we have done in this paper leads to 
a PdE on a graph with a large number of vertices. The grid points are the vertices of the graph, 
and the edges of the graph connect nearest neighbor grid points. Starting with an arbitrary graph, 
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FIGURE 10. A partial bifurcation diagram showing some of the solutions bifurcating 
from the 8-th and 10-th primary branches. Again, the dots at A = indicate 
solutions shown in Figures [13] and H31 The contour plots as a function of A are 
animated for the branches ending with the dots indicating symmetry types S\5 
(s7sl5.gif), S17 (s7sl7.gif), S 16 (s4sl6.gif), and S 2 2 (s4sl8s22.gif). 




Figure 11. A partial bifurcation diagram providing three additional symmetry 
types. For clarity, the trivial branch is not shown in this and the next figure. 
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FIGURE 12. A partial bifurcation diagram containing solutions of the seven remain- 
ing symmetry types. Primary branch 24 is the first branch with symmetry type S^- 
The symmetry types and Sn were found by searching the first one hundred 
primary branches, following only those branches which can lead to solutions with 
the desired symmetry. These two solutions are included for completeness, but their 
existence for the PDE would have to be confirmed with more modes and a higher 
level approximation of the eigenfunctions. 

our new suite of programs will analyze the symmetry of the graph and compute the bifurcation 
diagrams for the PdE on the graph. 

The programs we describe in the current paper will work with other superlinear odd / and other 
regions with hexagonal symmetry. The nonlinearity / needs to be superlinear since our program 
assumes that the branches eventually "go to the left." Our general program will not have this 
restriction; the GNGA and pmGNGA will be replaced by a single method of branch following that 
is able to go through fold points, and has no prejudice about the parameter increasing or decreasing. 
This new method of branch following has already been successfully implemented in [33]. We hope 
to write the new code so that a cluster of computers can be used in parallel, with each computer 
following a single branch at one time, under the control of a central PERL script. 

With minor modifications, our program would analyze the PDE ([I]) even when / is not odd. The 
appropriate bifurcation digraph for IS)q acting on L 2 (Q) is a subgraph of the digraph in Figure 
so the bifurcating branches would be followed properly unless the symmetry of the mother solution 
is incorrectly identified. The Perl scripts which start with the trivial branch would have to be 
modified, since u = is not a solution when / is not odd (unless /(0) = 0). If /(0) = 0, the 
trivial branch exists, but its bifurcations are not properly described by the bifurcation digraph in 
Figure [U and some special code would be needed to handle these bifurcations. 

It is valid to ask the question "does the GNGA converge" (as implemented in this current re- 
search). While we do not have a complete proof affirming the positive of this conjecture, many 
references contain relevant theorems. The GNGA is an implementation of Newton's method, which 
indeed converges under standard assumptions. In p3], one finds the classical fixed point itera- 
tion proof that Newton's method in l w converges when the initial guess is sufficiently close to a 
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Figure 13. The action of the generators of Bg on the plane, along with contour 
plots of solutions with symmetry types So, • • • , Sio at A = 0. Recall that S{ = [1^]. 
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FIGURE 14. Contour plots of solutions with symmetry types Su, . . . , S22 at A = 0. 
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Figure 15. A plot of «(2/27, 4^/27) function of the number of modes for 
the lowest energy solution at A = with symmetry type Sin. The point at M = 300 
matches the point labelled with Siq in Figure 



nondegenerate zero of the object function. This proof applies almost without change to the infi- 
nite dimensional case. Also addressed in [14] are algorithms where the object function and/or its 
derivative are only approximated; this would apply to our implementation due to numerical inte- 
gration errors, as well as owing to our imperfect knowledge of the eigenfunctions and corresponding 
eigenvalues. While not discussed exactly in the cited literature, elementary fixed point arguments 
indicate that the restriction of our object function V J to sufficiently large subspaces Bm will still 
result in convergent iterations. It would be worthwhile to string these type of results together 
in order to obtain a "best possible" GNGA convergence theorem. Monograph [13] gives an easy 
introduction into some of the details of implementing Newton's method to solve nonlinear prob- 
lems. Further, in the spirit of [7] and |35j . by the invariance of the Newton map, any convergence 
result should hold in fixed point subspaces corresponding to a given symmetry type. The articles 
[201 [35] and others by those authors discuss the convergence of algorithms similar to the GNGA, 
at times also considering symmetry restrictions. Finally, the well-known book \3>\ contains relevant 
convergence results for Newton and approximate Newton iterative fixed point algorithms. 

In summary, we have written a suite of programs that automatically computes the bifurcation 
diagram of the PDE (TjQ [2|). The program finds solutions with each of the 23 symmetry types by 
following solution branches which are connected to the trivial branch by a sequence of symmetry- 
breaking bifurcations. A thorough understanding of the possible symmetry-breaking bifurcations 
is required for this task. We introduced the bifurcation digraph, which summarizes the results of 
the necessary symmetry calculations. For the group Oq x Z2, these calculations were done by hand 
and verified by the GAP computer program [H [21] . In the future, we plan to implement automated 
branch following in systems where the symmetry group is so complicated that GAP is necessary. 
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